%% distShockInterventionMatrix.m
%
% Build Shock Intervention Matrices
%--------------------------------------------------------------------------

% Preliminaries
Mall = Nb*Na*Nz*Nr;
pointslistAll = (1:Mall)';

%1 = mortgage shock, 2 = liquidity shock, 3 = income shock
switch distShockType
case 1
    % Mortgage Shock
    assert(mortgageShock <= 0);
    a_shift = zeros(Na,1);
    a_wt_left = zeros(Na,1);
    for ai = 1:Na
        aprime = max(0,a(ai)+mortgageShock);
        aprime_left = discretize(aprime, a);   
        a_shift(ai) = ai - aprime_left;
        a_wt_left(ai) = 1 - ((aprime - a(aprime_left))./da);
    end
    %Now build end position for every point in grid
    %   Order is b,a,z,rm, and all that this should change is a
    ind_shift_left = a_shift.*Nb;
    ind_shift_left = reshape(transpose(repmat(ind_shift_left,1,Nb)),Na*Nb,1);
    ind_shift_left = repmat(ind_shift_left,Nz*Nr,1);
    ind_shift_right = ind_shift_left - Nb;

    a_wt_left_all = reshape(transpose(repmat(a_wt_left,1,Nb)),Na*Nb,1);
    a_wt_left_all = repmat(a_wt_left_all,Nz*Nr,1);

    SC = sparse(pointslistAll, pointslistAll-ind_shift_left, a_wt_left_all.*ones(Mall,1), Mall, Mall);
    SC = SC + sparse(pointslistAll, pointslistAll-ind_shift_right, (1-a_wt_left_all).*ones(Mall,1), Mall, Mall);    
    
    
case 2
    % Liquid Wealth Shock
    assert(liqShock > 0);
    b_shift = zeros(Nb,1);
    b_wt_left = zeros(Nb,1);
    for bi = 1:Nb
        bprime = min(bmax,b(bi)+liqShock);
        bprime_left = discretize(bprime, b);   
        b_shift(bi) = bprime_left - bi;
        b_wt_left(bi) = 1 - ((bprime - b(bprime_left))./db);
    end
    %Now build end position for every point in grid
    %   Order is b,a,z,rm, and all that this should change is b
    ind_shift_left = b_shift;
    ind_shift_left = repmat(ind_shift_left,Na*Nz*Nr,1);
    ind_shift_right = ind_shift_left + 1;

    b_wt_left_all = repmat(b_wt_left,Na*Nz*Nr,1);

    SC = sparse(pointslistAll, pointslistAll+ind_shift_left, b_wt_left_all.*ones(Mall,1), Mall, Mall);
    SC = SC + sparse(pointslistAll, pointslistAll+ind_shift_right, (1-b_wt_left_all).*ones(Mall,1), Mall, Mall);
    
    
case 3
    % Income Shock (High -> Medium; Medium -> Low) 
    assert(incomeRelativeWeightShiftL > 0);
    assert(incomeRelativeWeightShiftM > 0);
    assert(incomeRelativeWeightShiftH > 0);
    
    SC = sparse(Mall,Mall);
    for rind = 1:Nr
        SCsmall = sparse(Nb*Na*Nz,Nb*Na*Nz); 
        %Low income stays put
        SCsmall = SCsmall + sparse([1:Nb*Na], [1:Nb*Na], ones(1, Nb*Na), Nb*Na*Nz,Nb*Na*Nz);
        %Adjust middle income state to low income state
        SCsmall = SCsmall + sparse([Nb*Na+1:Nb*Na*2],[Nb*Na+1:Nb*Na*2], (1-incomeRelativeWeightShiftM)*ones(1,Nb*Na), Nb*Na*Nz,Nb*Na*Nz);
        SCsmall = SCsmall + sparse([Nb*Na+1:Nb*Na*2],[1:Nb*Na], incomeRelativeWeightShiftM*ones(1,Nb*Na), Nb*Na*Nz,Nb*Na*Nz);
        %Adjust high income state to middle income state
        SCsmall = SCsmall + sparse([Nb*Na*2+1:Nb*Na*3],[Nb*Na*2+1:Nb*Na*3], (1-incomeRelativeWeightShiftH)*ones(1,Nb*Na), Nb*Na*Nz,Nb*Na*Nz);
        SCsmall = SCsmall + sparse([Nb*Na*2+1:Nb*Na*3],[Nb*Na+1:Nb*Na*2], incomeRelativeWeightShiftH*ones(1,Nb*Na), Nb*Na*Nz,Nb*Na*Nz);

        %Combine across interest rates
        SC(Nb*Na*Nz*(rind-1)+1:Nb*Na*Nz*rind,Nb*Na*Nz*(rind-1)+1:Nb*Na*Nz*rind)=SCsmall;
    end
    clearvars SCsmall;
end

